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A new algorithm for simulation of theories with dynamical fermions is presented. The algorithm is based on 
obtaining the new configuration U from the old one U by solving the equation M{U)rj = uM{U)ri, where M is 
fermionic operator, r/ is random Gaussian vector, and ui is random real number close to unity. This algorithm can 
be used for acceleration of current simulations in theories with fermions. 



From the earliest days of Monte Carlo simula- 
tions fermionic fields have caused annoying diffi- 
culty, which stems from their being anticommut- 
ing variables. Although for the most actions in 
use one can eliminate fermions by an analytic in- 
tegration, the resulting expressions involve deter- 
minants of large matrices, making the numerical 
simulations expensive. 

Over the years many tricks have been devel- 
oped to solve this problem. Hybrid Monte Carlo 
(HMC) is now considered to be the standard 
simulation tool. In recent years it was challenged 
by the Multiboson (MB) method. The most pop- 
ular version of MB, proposed by M. Liischer 
has been intensively studied by many authors (see 
|4| and references therein). The MB algorithm, 
proposed by A. Slavnov |^ , is less known, and was 
tested in Besides HMC and MB, let me men- 
tion some interesting ideas, based on the poly- 
mer [Q and Jordan- Wigner occupation number 
|8| representations, the direct evaluation of Grass- 
mann integrals |^ , the separation of low and high 
eigenmodes of Dirac operator jl^ , and the direct 
simulation of loop expansion by means of using 
the stochastic estimators |]lT| . 

Despite the variety of different approaches to 
the problem, there is a common impression that 
all existing algorithms remain inefficient [p^ . 
This is clearly seen if we compare the compu- 
tational cost for the theories with dynamical 
fermions from the one side, and purely bosonic 
theories from the other side. Such a situation 
suggests that one should not stop the attempts 



to obtain relatively cheap fermion algorithm. 

In this contribution I present a new computa- 
tional strategy for treating dynamical fermions in 
Monte Carlo simulations. It has the virtues of ex- 
actness, nonlocality and finite step-size of update 
Unlike the other exact algorithms, the CPU-time 
required for the update in this algorithm grows 
only linearly with the volume V of the system. 

Suppose that we aim at sampling the partition 
function of theory with two flavors of degenerate 
fermion fields: 

Z,.„,^lDe,{M<mm)^U (1) 

where M is fermion operator, U denotes the 
bosonic fields coupled to fermions. (For simplic- 
ity of formulae I do not include the purely bosonic 
contribution.) Then starting from old configura- 
tion U, one executes the following instructions: 

• generate random vector i] with Gaussian 
distribution: PgM = Z^^-^^^^ 



• generate uj £ 
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• find the new configuration U by solving the 
equation: 



Mr] — LoAIr] 



(2) 



where < e < 1 is an algorithmic parameter. 
Here and below I use the following short nota- 
tions: M = M[U]: M = M[U] . 
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In Ref. |l[ it was proved that as far as the so- 
lution of eq. (H) is found in a reversible way Q, the 
algorithm satisfies the detailed balance condition. 

It is clear that choosing the procedure of solv- 
ing the eq.(||) is crucial in this approach. Firstly, 
because obtaining the solution U is the only 
computationally expensive ingredient of the al- 
gorithm. Secondly, because not every procedure 
of solving the eq. ^ can provide ergodicity. Find- 
ing a good procedure of solving the eq. (||) presents 
a serious challenge for the future investigations. 
The potential gain worth the efforts: e.g. finding 
some analytic solutions can give us a very cheap 
fermionic algorithm, comparable in cost with the 
algorithms for purely bosonic theories. 

Making the first attempt to test the new algo- 
rithm, I solved the eq.(||) numerically by using 
the local iterative minimization of the quantity 

R=\{M -LuM)rj\'^ (3) 

for fixed U,Lu,r]. Starting from U = U, the mini- 
mization proceeds until R < 5 in reached, where 

5 determines the accuracy of solution. 

The tests were done for SU(3) QCD at 8'* lat- 
tice for the partition function (|l|) with 

M[U] - 1 - PD^oDoe (4) 

where D[U] is Wilson difference operator. The 
even-odd preconditioning was used to reduce the 
computational cost. The hopping parameter was 
k = 0.2, which gave the plaquette V = 0.0089(1). 
The minimization of functional (H) was imple- 
mented by making the random moves in each 
color direction for all links lexicographically. 

I checked the existence of solutions of eq.(||) 
at w £ [0.9, 1.1] for typical equilibrated config- 
urations U. It was found that for considered lo 
one always finds some solution. On the aver- 
age solving the eq.(||) at e = 0.05 with precision 

6 — 5 * 10~^ required 42 minimization iterations 
(it was checked, that improving the precision did 
not affect the results) . Since the cost of one iter- 
ation is roughly equal to 2 * Ngenerator = 16 mul- 
tiplications by matrix M, the total cost of finding 

-"^i.e. the probability to obtain configuration U starting 
from U at any ui should be equal to probability to obtain 
U starting from U at 1/u! 



the solution was approximately 670 matrix mul- 
tiplication. This has to be compared with the 
cost for generating one trajectory with HMC. At 
step-size r = 0.033 and trajectory length equal 
to 1 (this ensures 70% acceptance), one trajec- 
tory cost approximately 4420 matrix multiplica- 
tions, i.e. ~ 6.6 times more expensive than solv- 
ing eq.(|). 

My conclusion was that the usage of minimiza- 
tion of functional (||) for solving the eq.(||) is not 
favorable, because the algorithm became noner- 
godic. The problem is that probability V[uj,r]] 
to accept w < 1 is strongly suppressed by large 
factor in the exponent, if uj is far from 1. In- 
deed, the squared norm of Gaussian vector can 
be estimated as ~ Idof, where Idof is the 

dimensionality of vector space on which i] is de- 
fined. For even-odd preconditioned SU(3) QCD 
at 8^ lattice one has Idof = 24576! Therefore, 
using the new algorithm alone, one is restricted 
to choose between two unfavorable possibilities: 
1) Using e ^ 1/Idof- It was observed that at 
such small values of e the evolution of U fields be- 
comes very slow, because the new configuration 
U always lies too close to the old one; 2) Using 
e 3> 1/Idof- Then w < 1 is almost never accepted, 
and the configurations are becoming smother. 

In the second case the trouble appears due to 
the nonergodicity of algorithm at large e. Never- 
theless, even at large e our algorithm can be used 
for acceleration of other fermionic algorithms like 
HMC or MB, since the ergodicity is provided by 
them. 

I tested the combination of the new algorithm 
and HMC. After each HMC trajectory the global 
move in configuration space was implemented by 
using our algorithm at e = 0.05. The efficiency of 
this algorithmic mixture was estimated by mea- 
suring the autocorrelation times for the plaque- 
tte. The resuhs were: 1) Pure HMC, 1000 tra- 
jectories: nnt = 1-4(2); 2) HMC-fnew alg., 1000 
(traj -I- w-update): Tint = 0.7(1). Note, that au- 
tocorellation time was reduced almost at no cost, 
because the global update of the new algorithm is 
much cheaper, than HMC trajectory. Of course, 
these tests are not very illustrative, since the au- 
tocorrelation times for theory with pure fermionic 
term at this kappa value are rather small. The 
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further tests in models of practical relevance are 
needed. 

Let me emphasize, that nonergodicity of the 
new algorithm was produced by the particular 
way of solving the eq.(||), which forces us to use 
large e. One should think of finding some other 
procedure, which allows to use e = (or e « 0). 
Unfortunately, in this case the search for the so- 
lution of eq.(^ by minimizing the functional (|^) 
does not work, because one immediately gets the 
trivial solution U — U. However, in typical case 
(e.g. SU(3) QCD), the solutions of eq.(^ are 
highly degenerate |^ and the subspace of non- 
trivial solutions is not empty. Finding the re- 
versible procedure of nontrivial solving the equa- 
tion Mr] = Mr], which samples the configuration 
space fast enough, is the perspective subject for 
future investigations. 

I conclude that the efficiency of the new fermion 
algorithm drastically depends on the clever choice 
of procedure of solving the eq.(^). Finding the 
new efficient way of solving this equation can pro- 
vide a very cheap simulation method. 
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